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The nature of the unknown sources of ultra-high energy cosmic rays can be revealed 
through the detection of the GZK feature in the cosmic ray spectrum, resulting from 
the production of pions by ultra-high energy protons scattering off the cosmic mi- 
crowave background. Here we show that the GZK feature cannot be accurately 
determined with the small sample of events with energies ~ 10 20 eV detected thus 
far by the largest two experiments, AGASA and HiRes. With the help of numer- 



ical simulations for the propagation of cosmic rays, we find the error bars around 



the GZK feature are dominated by fluctuations which leave a determination of the 
GZK feature unattainable at present. In addition, differing results from AGASA and 
HiRes suggest the presence of ~ 30% systematic errors that may be due to discrep- 
ancies in the relative energy determination of the two experiments. Correcting for 
these systematics, the two experiments are brought into agreement at energies be- 
low ~ 10 20 eV. After simulating the GZK feature for many realizations and different 
injection spectra, we determine the best fit injection spectrum required to explain 
the observed spectra at energies above 10 18 ' 5 eV. We show that the discrepancy 
between the two experiments at the highest energies has low statistical significance 
(at the 2 a level) and that the corrected spectra are best fit by an injection spectrum 
with spectral index ~ 2.6. Our results clearly show the need for much larger exper- 
iments such as Auger, EUSO, and OWL, that can increase the number of detected 
events by 2 orders of magnitude. Only large statistics experiments can finally prove 
or disprove the existence of the GZK feature in the cosmic ray spectrum. 
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1 Introduction 



The presence or lack of a feature in the spectrum of ultra-high energy cos- 
mic rays (UHECRs) is key in determining the nature of these particles and 
their sources. Astrophysical proton sources distributed homogeneously in the 
universe produce a feature in the spectrum due to the production of pions off 
the cosmic microwave background. This feature, consisting of a rather sharp 
suppression of the flux, occurs at energies above 7 x 10 19 eV, as a result of 
the threshold in the production of pions in the final state of a proton-photon 
inelastic interaction. This important prediction was made independently by 
Greisen and by Zatsepin and Kuzmin [1]. The resulting spectral feature is 
now known as the GZK cutoff (or feature, as we prefer to call it). Alternative 
models for UHECR sources that involve new physical processes may evade 
the presence of this feature (see, e.g., [2]). Recent reviews on the origin and 
propagation of the ultra-high energy cosmic rays can be found in [3,4], while 
a recent review of the observations can be found in [5]. 

The detection of cosmic ray events with energy above E GZK ~ 7x 10 19 eV does 
not necessarily imply that the GZK feature is not present: what characterizes 
the presence of the GZK feature is the relative number of events above and be- 
low Eqzk when both sides of the spectrum can be accurately determined. The 
steep injection spectra required to fit the observations below Eqzk imply that 
only a handfull of events above 10 20 eV can be detected during the operation 
time of experiments such as AGASA and HiRes. This makes the identification 
of the GZK feature by these experiments extremely difficult. The problem 
is exacerbated by the fluctuations due to the discreteness of the process of 
photo-pion production, as will be discussed below. These uncertainties need 
to be considered when attempting a determination of the best fit injection 
spectrum of the particles, and in order to quantify the statistical significance 
of the presence or absence of the GZK feature in the observed spectrum. 

The discrepancy between the results of the two largest experiments has gener- 
ated much debate. AGASA [6] reports a higher number of events above E GZ k 
than expected while HiRes[7-9] reports a flux consistent with the GZK feature. 
Here we investigate in detail the statistical significance of this discrepancy as 
well as the significance of the presence or absence of the GZK feature in the 
data. We find that neither experiment has the necesseray statistics to estab- 
lish if the spectrum of UHECRs has a GZK feature. In addition, a systematic 
error in the energy determination of the two experiments seems to be required 
in order to make the two sets of observations compatible in the low energy 
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range, 10 18 ' 5 — 10 196 eV, where enough events have been detected to make the 
measurements reliable. The combined systematic errors in the energy deter- 
mination is likely to be ~ 30%. If we decrease the AGASA energies by 15% 
while increasing HiRes energies also by 15%, the two experiments predict com- 
patible fluxes at energies below Eqzk and at energies above Eqzk the fluxes 
are within ~ 2a of each other. In this case, the best fit injection spectrum 
has a spectral index of ~ 2.6 but a determination of the GZK feature has 
very low significance. The detection or non-detection of the GZK feature in 
the cosmic ray spectrum remains open to investigation by future generation 
experiments, such as the Pierre Auger project [10] and the EUSO [11] and 
OWL [12] experiments. 

This paper is planned as follows: in §2 we describe our simulations. In §3 we 
illustrate the present observational situation, limiting ourselves to AGASA 
and HiResI, and compare the data to the predictions of our simulations. We 
conclude in §4. 



2 UHECR spectrum simulations 

We assume that ultra-high energy cosmic rays are protons injected with a 
power-law spectrum in extragalactic sources. The injection spectrum is taken 
to be of the form 



where 7 is the spectral index, a is a normalization constant, and E max is the 
maximum energy at the source. Here we fix E max = 10 215 eV, large enough not 
to affect the statistics at much lower energies. As shown in [13], the induced 
spectrum of a uniform distribution of sources in space is almost indistinguish- 
able from a distribution with the observed large scale structure in the galaxy 
distribution. Based on this result, we assume a spatially uniform distribution 
of sources and do not take into account luminosity evolution in order to avoid 
the introduction of additional parameters. 

We simulate the propagation of protons from source to observer by including 
the photo-pion production, pair production, and adiabatic energy losses due to 
the expansion of the universe. In each step of the simulation, we calculate the 
pair production losses using the continuous energy loss approximation given 
the small inelasticity in pair production (2m e /m p ~ 10~ 3 ). For the rate of 
energy loss due to pair production at redshift z = 0, (3 PP (E, z = 0), we use the 
results from [14,15]. At a given reshift z > 0, 



F(E)dE = aE-i expi-E/EjnJdE 



(1) 



(3 PP (E, z) = (l + z) 3 (3 pp ((l + z)E, z = 0). 



(2) 
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Similarly, the rate of adiabatic energy losses due to redshift is calculated in 
each step using 



PME, z) 



n M (i + z) 



1/2 



(3) 



with H = 75 km s x Mpc \ 

The photo-pion production is simulated in a way similar to that described in 
ref. [13]. In each step, we first calculate the average number of photons able 
to interact via photo-pion production through the expression: 



where l(E, z) is the interaction length for photo-pion production of a proton 
with energy E at redshift z and As is a step size, chosen to be much smaller 
than the interaction length (typically we choose As = 100 kpc/(l + z) 3 ). 

In Fig. 1 we plot the interaction length for photopion production used in [3] 
(solid thin line), and in [16] (triangles). The dashed line is the result of our 
calculations (see below), which is in perfect agreement with the results of 
[3,16]. The apparent discrepancy at energies below 10 19 5 eV with the predic- 
tion of Ref. [3] is only due to the fact that we consider only microwave photons 
as background, while in [3] the infrared background was also considered. For 
our purposes, this difference is irrelevant as can be seen from the loss lengths 
plotted in Fig. 1. The rightmost thick solid line is the loss length for photo- 
pion production [3], while the other thick solid line is the loss length for pair 
production. In the present calculations, we do not use the loss length of pho- 
topion production which is related to the interaction length through an angle 
averaged inelasticity. Unlike what was done in [13], in the current simulations 
we evaluate the inelasticity for each single proton-photon scattering using the 
kinematics, rather than adopting an angle averaged value. 

We calculate the interaction length, 1(E), as: 

1(E)' 1 = fden(e) + [d^-^a(s) (5) 



where n(e) is the number density of the CMB photons per unit energy at 
energy e, [3 is the velocity of the proton, /j, is the cosine of the interaction 
angle, and a(s) is the total cross section for photo-pion production for the 
squared center of mass energy 

s = m 2 + 2eE (1 - fip) . (6) 
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For the calculation of the interaction length we adopt the p + 7 — > hadrons 
cross section given in [17], considering only the photons of the cosmic mi- 
crowave background as targets. The calculated interaction length (see dashed 
line in Fig. 1) is in good agreement with the interaction length calculated in 
[3] and in [16]. 

Once the interaction length is known, we then sample a Poisson distribution 
with mean (N p h(E, As)), to determine the actual number of photons encoun- 
tered during the step As. When a photo-pion interaction occurs, the energy 
e of the photon is extracted from the Planck distribution, n ph (e,T(z)), with 
temperature T(z) —T (l + z), where T = 2.728 K is the temperature of the 
cosmic microwave background at present. Since the microwave photons are 
isotropically distributed, the interaction angle, 9, between the proton and the 
photon is sampled randomly from a distribution which is flat in /i — cos9. 
Clearly only the values of e and 9 that generate a center of mass energy above 
the threshold for pion production are considered. The energy of the proton in 
the final state is calculated at each interaction from kinematics. The simula- 
tion is carried out until the statistics of events detected above some energy 
reproduces the experimental numbers. By normalizing the simulated flux by 
the number of events above an energy where experiments have high statis- 
tics, we can then ask what are the fluctuations in numbers of events above 
a higher energy where experimental results are sparse. The fits are therefore 
most sensitive to the energy regions below E GZ k and give a good estimate of 
the uncertainties in the present experiments for energies above Eqzk- In this 
way we have a direct handle on the fluctuations that can be expected in the 
observed flux due to the stochastic nature of photo-pion production and to 
cosmic variance. 

The simulation proceeds in the following way: a source distance is generated at 
random from a uniform distribution in a universe with = 0.7 and Q m = 0.3. 
In a Euclidean universe, the flux from a source would scale as r~ 2 where r is the 
distance between the source and the observer. On the other hand, the number 
of sources between r and r + dr would scale as r 2 , so that the probability that 
a given event has been generated by a source at distance r is independent of 
r: sources at different distances have the same probability of generating any 
given event. In a flat universe with a cosmological constant, this is still true 
provided the distance r is taken to be 



where t g is the age of the universe when the event was generated, to is the 
present age of the universe, and R(t) is the scale factor of the universe. Once a 
source distance has been selected at random, a particle energy is also assigned 




(7) 
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Fig. 1. Interaction length for photopion production as calculated in this paper 
(dashed line) compared to the interaction length of [3] (solid thin line) and of [16] 
(triangles). The thick solid lines are the loss lengths for photopion production (on 
the right) and of pair production (on the left). 

from a distribution that reflects the injection spectrum, chosen as in Eq. (1). 
This particle is then propagated to the observer and its energy recorded. This 
procedure is repeated until the number of events above a threshold energy, E th 
is reproduced. With this procedure we can assess the significance of results 
from present experiments with limited statistics of events. There is an addi- 
tional complication in that the aperture of the experiment usually depends 
on energy. This is taken into account by allowing the event to be detected 
or not detected depending upon the function H(E) that describes the energy 
dependence of the aperture. 

We only study the spectrum above 10 18 ' 5 eV where the flux is supposed to be 
dominated by extragalactic sources. For this energy range, we focus on the ex- 
periments that have the best statistics: AGASA and HiResI. For the AGASA 
experiment the exposure is basically flat above 10 19 eV, while for HiRes the ex- 
posure is as plotted in Fig. (2) [9]. For AGASA data, the simulation is stopped 
when the number of events above E t h = 10 19 eV equals 866. For HiRes this 
number is 300. Note that while for AGASA the number of detected events ac- 
tually corresponds to the generated events, for HiRes the number of detected 
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Fig. 2. Aperture of the HiResI experiment as a function of the energy from [9]. 

events requires a correction due to the energy dependence of the aperture 
H(E). This correction allows one to reconstruct the observed spectrum. The 
statistical error in the energy determination is accounted for in our simulation 
by generating a detection energy chosen at random from a Gaussian distribu- 
tion centered at the arrival energy E and with width AE/E = 30% for both 
experiments. 

Our simulations reproduce well the predictions of analytical calculations, in 
particular at the energies where energy losses may be approximated as con- 
tinuous. In Fig. 3, we compare the results of our simulation with analytical 
calculations carried out as in ref. [18]. The curves plotted in the figure are the 
so called modification factors, defined in [18,19] for three different values of 
the source redshift (z = 0.002, z = 0.02 and z = 0.2 from top to bottom). The 
differential injection spectrum is taken to be a power law E~ 21 . The points 
with error bars are the results of our simulations with 2 x 10 6 particles pro- 
duced by sources at the redshifts given above. The agreement between the 
simulated and analytical calculations at low energies is clear. At the ener- 
gies where photo-pion production becomes important, simulations predict a 
slightly larger flux than analytical calculations. This is a well known result, 
and is due to the discreteness of photo-pion energy losses, that allow some 
particles to reach the detector without appreciable losses. 

In Fig. 4 we compare the results of our simulations (points with error bars) 
with analytical calculations of the diffuse flux of UHECRs (continuous line) 
expected if sources with no luminosity evolution inject UHECRs with a spec- 
trum E~ 2 - 7 [20,21]. The agreement is evident. 
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log 10(E) [eV] 



Fig. 3. Comparison between analytical calculations and the results of our simulation 
for the modification factor, for injection spectrum E~ 21 [18] and three values of the 
source redshift (z = 0.002, z = 0.02 and z = 0.2 from top to bottom). 
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Fig. 4. Comparison between the results of our simulation and the analytical calcu- 
lations of [20,21] for the case of astrophysical sources injecting cosmic rays with a 
spectrum E~ 2 ' 7 and no luminosity evolution. 

3 AGASA versus HiResI 



The two largest experiments that measured the flux of UHECRs report ap- 
parently conflicting results. The data of AGASA and HiResI on the flux of 
UHECRs multiplied as usual by the third power of the energy are plotted in 
Fig. 5 (the squares are the HiResI data while the circles are the AGASA data). 
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20.5 
loglO(E) [eV] 



Fig. 5. Circles show the AGASA spectrum from [6] while squares show the HiResI 
spectrum from [9]. 

Table 1 

Number of events for AGASA and HiResI detected above the energy thresholds 
reported in the first column. 



log(E th ) (eV) 


AGASA 


HiResI 




data 


-15% 


data +15% 


19 


866 


651 


300 367 


19.6 


72 


48 


27 39 


20 


11 


7 


1 2.2 



The figure shows that HiResI data are systematically below AGASA data and 
that HiResI sees a suppression at ~ 10 20 eV that resembles the GZK feature 
while AGASA does not. 



We apply our simulations here to the statistics of events of both AGASA 
and HiResI in order to understand whether the discrepancy is statistically 
significant and whether the GZK feature has indeed been detected in the 
cosmic ray spectrum. The number of events above 10 19 eV, 10 196 eV and 10 20 
eV for AGASA and HiResI are shown in Table 1. For reasons that will be 
clear below, we also show in the same table the number of events above the 
same energy thresholds, in the case that AGASA has a systematic error that 
overestimates the energy by 15% while HiResI systematically underestimates 
the energy by 15%. In order to understand the difference, if any, between 
AGASA and HiRes data we first determine the injection spectrum required 
to best fit the observations. In order to do this, we run 400 realizations of the 
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Table 2 

X 2 for fits to AGASA and HiResI data above 10 18 - 5 eV, 10 19 eV, and 10 19 - 6 eV for 
varying spectral index 7. In parenthesis the number of degrees of freedom. 







AGASA 






HiResI 




7 


Xi8 5(!7) 


Xiq(12) 


X?q fi(6) 


Xi8^15) 

'VlO.O \ / 


Xw(10) 

/v. 1 y v / 


/viy.o \ / 


2.3 


1694 


34 


17 


145 


29 


23 


2.4 


1215 


24 


12 


80 


19 


15 


2.5 


724 


19 


10 


36 


14 


11 


2.6 


248 


16 


10 


14 


9 


7 


2.7 


82 


17 


11 


33 


7 


5 


2.8 


80 


21 


13 


126 


7 


4 


2.9 


316 


27 


15 


257 


9 


4 



AGASA and HiRes event statistics, as reported in the column labelled data 
in Table 1. 

The injection spectrum is taken to be a power law with index 7 between 2.3 and 
2.9 with steps of 0.1. For each injection spectrum we calculated the x 2 indicator 
(averaged over 400 realizations for each injection spectrum). The errors used 
for the evaluation of the x 2 are due to the square roots of the number of 
observed events. As reported in Table 2, the best fit injection spectrum can 
change appreciably depending on the minimum energy above which the fit is 
calculated. In these tables we give Xg(iV), where N is the number of degrees 
of freedom, while the subscript, e, is the logarithm of E t h (in eV), the energy 
above which the fit has been calculated. The numbers in bold-face correspond 
to the best fit injection spectrum. These fits are dominated by the low energy 
data rather than by the poorer statistics at the higher energies. 

If the data at energies above 10 18 ' 5 eV are taken into account for both ex- 
periments, the best fit spectra are E~ 2 ' 8 for AGASA and E~ 2S for HiRes. If 
the data at energies above 10 19 eV are used for the fit, the best fit injection 
spectrum is E~ 2 - 6 for AGASA and between E~ 2 - 7 and E~ 2 - 8 for HiRes. If the 
fit is carried out on the highest energy data (E > 10 196 eV), AGASA prefers 
an injection spectrum between E~ 2,b and E~ 2S , while E~ 2,8 or E~ 2M fit better 
the HiRes data in the same energy region. Note that the two sets of data un- 
corrected for any possible systematic errors require different injection spectra 
that change with E th . 

In order to quantify the significance of the detection or lack of the GZK 
flux suppression, we report in Tables 3 and 4 the mean number of events 
above the indicated energy threshold, (N(E > E t h)), for different injection 
specta. In parenthesis , we show the discrepancy between the observations as 
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Table 3 

Number of events expected above E t h(eV) for different injection spectra assuming 
the AGASA statistics above 10 19 eV. In parenthesis are the number of standard de- 
viations, a, between the observed and expected number of events. In square brackets 
are the discrepancies calculated with a combined error bar of simulation and obser- 
vation uncertainties, atot- 





7 = 2.5 


7 = 2.6 


7 = 2.8 


1Q 19.6 


65 ±8.2 (+0.5) [+0.3] 


58 + 7.6 (+1.4) [+1.0] 


46 ±6.8 (+2.8) [+2.2] 


10 20 


3.5 ± 1.9 (+2.4) [+2.1] 


2.8 + 1.7 (+2.6) [+2.3] 


2.0 + 1.4 (+2.8) [+2.6] 



Table 4 

Number of events expected above E t h(eV) for different injection spectra assuming 
the HiResI statistics above 10 19 eV from Table 1. In parenthesis are the number of 
a between the observed and expected number of events. In square brackets are the 
number of atot ■ 



Eth 


7 = 2.6 


7 = 2.7 


7 = 2.8 


1Q 19.6 


31 ±5.6 (-0.8) [-0.6] 


28 + 5.3 (-0.2) [-0.1] 


26 ±5.2 (+0.3) [+0.2] 


1Q 20 


1.9 + 1.4 (-0.9) [-0.5] 


1.5 + 1.2 (-0.5) [-0.3] 


1.3 + 1.2 (-0.3) [-0.2] 



in Table 1 and the simulations in number of standard deviations, a, where a 2 = 
(N(E > E th f - (N(E > E th )) 2 ). From Tables 3 and 4 we can see that while 
HiResI is consistent with the existence of the GZK feature in the spectrum of 
UHECRs, AGASA detects an increase in flux, but only at the ~ 2.5a level for 
the best fit injection spectra. 



A more graphical representation of the uncertainties involved are displayed in 
Figs. 6 for AGASA and 7 for HiRes, for two choices of the injection spectrum. 
These plots show clearly the low level of significance that the detections above 
Eqzk have with low statistics. The large error bars that are generated by 
our simulations at the high energy end of the spectrum are mainly due to the 
stochastic nature of the process of photo-pion production: in some realizations 
some energy bins are populated by a few events, while in others the few par- 
ticles in the same energy bin do not produce a pion and get to the observer 
unaffected. The large fluctuations are unavoidable with the extremely small 
statistics available with present experiments. On the other hand, the error 
bars at lower energies are minuscule, so that the two data sets (AGASA and 
HiResI) cannot be considered to be two different realizations of the same phe- 
nomenon. Instead, systematic errors in at least one if not both experiments 
are needed to explain the discrepancies at lower energies. 

Taking into account the (theoretical) error bars in the analysis makes the 
significance of the presence or absence of the GZK feature much weaker than 
the consistency checks shown in Tables 3 and 4. In order to estimate this effect, 
we proceed in the following way: we calculate the expected number of events 
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AGASA 


y=2.6 


E>10" 


866 


E>10" 6 


57.49 ± 7.59 


E>10 2 ° 


2.82 ± 1.71 



20.5 21 
loglO(E) [eV] 



AGASA 


7=2.8 


E>10" 


866 


E>10 1 "' 


45.76 ± 6.77 


E>10 20 


2.02 ± 1.35 




20.5 21 
log 10(E) [eV] 



Fig. 6. Simulations of AGASA statistics with injection spectra E~ 2,e (upper plot) 
and E~ 2,8 (lower plot). The crosses with error bars are the results of simulations, 
while the grey points are the AGASA data. 

above some threshold with its corresponding standard deviation (<J S i m ), as 
determined by the fluctuations in the flux simulation. The observed number 
of events above the same threshold is also known with the error bar a b s - 
The discrepancy between the two is now calculated using the error a to t — 
( a sim + a obsY^ 2 - O ur resUi ts are summarized in Tables 3 for AGASA and 
4 for HiRes. The numbers with error bars are the simulated expectations, 
while the discrepancy between simulations and observations, calculated as 
described above is in square brackets, in units of atot- If becomes clear that 
the effective discrepancy between predictions and the AGASA data is at the 
level of 2.1 — 2.5a. Therefore a definite answer to the question of whether the 
GZK feature is there or not awaits a significant improvement in statistics at 
high energies. 

As seen in Fig. 5, the difference between the AGASA and HiResI spectra is 
not only in the presence or absence of the GZK feature: the two spectra, when 
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E>10 15 
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: 300 
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: 1.50 


k 1.23 



18.5 



19 



19.5 



20 



20.5 21 
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Fig. 7. Simulations of HiRes statistics with injection spectra E 2,6 (upper plot) and 
E~ 2 - 7 (lower plot). The crosses with error bars are the results of simulations, while 
the squares are the HiRes data. 

multiplied by E 3 , are systematically shifted by about a factor of two. This 
shift suggests that there may be a systematic error either in the energy or the 
flux determination of at least one of the two experiments. Possible systematic 
effects have been discussed in [22] for the AGASA collaboration and in [9] for 
HiResI. At the lower end of the energy range the offset may be due to the 
notoriously difficult determination of the AGASA aperture at threshold while 
the discrepancies at the highest energies is unclear at the moment. In any 
case, a systematic error of ~ 15% in the energy determination is well within 
the limits that are allowed by the analysis of systematic errors carried out by 
both collaborations. 



In order to illustrate the difficulty in determining the existence of the GZK 
feature in the observed data in the presence of systematic errors, we split 
the energy gap by assuming that the energies as determined by the AGASA 
collaboration are overestimated by 15%, while the HiRes energies are underes- 
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Table 5 

X 2 for AGASA and HiResI in which a correction for a systematic 15% overestimate 
of the energies has been assumed for AGASA and a 15% underestimate of the 
energies has been assumed for HiResI. 







AGASA 






HiResI 




7 


Xke(15) 


X? 9 (H) 


X?9.6(5) 


X? 8 . 6 (14) 


x? 9 (io) 


X? 9 . 6 (4) 


2.3 


505 


18 


12 


79 


13 


7 


2.4 


351 


13 


8.5 


40 


7 


4 


2.5 


188 


9 


5.6 


13 


3.7 


2.0 


2.6 


54 


9 


5.6 


6 


2.0 


1.1 


2.7 


20 


11 


6.4 


23 


3.1 


1.4 


2.8 


54 


15 


7.2 


94 


6 


2.4 


2.9 


145 


20 


9.1 


176 


10 


4 



timated by the same factor. The number of events above an energy threshold 
is again reported in Table 1, in the column labelled 15%. In this case the total 
number of events above 10 19 eV is reduced for AGASA from 866 to 651, while 
for HiResI it is enhanced from 300 to 367. We ran our simulations with these 
new numbers of events and repeat the statistical analysis described above. 
The values of x 2 obtained in this case are reported in Table 5. 

For AGASA, the best fit injection spectrum is now between E~ 2 - 5 and E~ 2 - 6 
above 10 19 eV and above 10 19 6 eV (the x 2 values are identical). For the HiRes 
data, the best fit injection spectrum is E~ 2S for the whole set of data, inde- 
pendent of the threshold. It is interesting to note that the best fit injection 
spectrum appears much more stable after the correction of the 15% system- 
atics has been carried out. Moreover, the best fit injection spectra as derived 
for each experiment independently coincides for the corrected data unlike the 
uncorrected case. This suggests that combined systematic errors in the energy 
determination at the ~ 30% level may in fact be present. 

The expected numbers of events with energy above 10 196 eV and above 10 20 
eV with the deviation from the data are reported in Tables 6 and 7: while 
HiResI data remain in agreement with the prediction of a GZK feature, the 
AGASA data seem to depart from such prediction but only at the level of 
~ 1.8(7. The systematics in the energy determination further decreased the 
significance of the GZK feature, such that the AGASA and HiResI data are 
in fact only less than 2a away from each other. 

We can use the same procedure illustrated above to evaluate the effect of the 
error bars in the simulated data compared to the data corrected by 15%. The 
results are reported in square brackets in Tables 6 (for AGASA) and 7 (for 
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Table 6 

Number of events expected above E t h (eV) for AGASA energies decreased sys- 
tematically by 15%. In parenthesis is the number of standard deviations between 
observations and simulations, a. In square brackets are the discrepancies calculated 
in units of a tot- 



Eth 


7 = 2.5 


7 = 2.6 


7 = 2.7 


10 19.6 


49 ±6.9 (+0.2) [+0.1] 


43 ±6.5 (+0.8) [+0.5] 


39 ±6.1 (+1.3) [+1.0] 


1Q 20 


2.6 + 1.6 (+1.7) [+1.4] 


2.3 + 1.5 (+1.8) [+1.5] 


1.8 + 1.4 (+2.0) [+1.7] 



Table 7 

Number of events expected above E t h (eV) for HiResI energies increased systemat- 
ically by 15%. In parenthesis is the number of standard deviations between obser- 
vations and simulations, a. In square brackets are the discrepancies calculated in 
units of a tot- 



Eth 


7 = 2.5 


7 = 2.6 


10 19.6 


43 ±6.3 (-0.6) [-0.4] 


38 ±6.0 (+0.1) [+0.1] 


10 20 


2.8 + 1.7 (-0.4) [-0.3] 


2.3 + 1.5 (-0.1) [-0.1] 



HiRes), showing that the effective discrepancy between expectations (with 
uncertainties due to discrete energy losses and cosmic variance) and AGASA 
data is even smaller, only at the 1.5<r level. In Fig. 8, we plot the simulated 
spectra for injection spectrum E~ 2M and compare them to observations of 
AGASA (upper plot) and HiRes (lower plot). 



4 Conclusions 

We considered the statistical significance of the UHECR spectra measured by 
the two largest experiments in operation, AGASA and HiRes. The spectrum 
released by the HiResI collaboration seems to suggest the presence of a GZK 
feature. This has generated claims that the GZK cutoff has been detected, 
reinforced by data from older experiments [23]. However, no evidence for such 
a feature has been found in the AGASA experiment. We compared the data 
with theoretical predictions for the propagation of UHECRs on cosmological 
distances with the help of numerical simulations. We find that the very low 
statistics of the presently available data hinders any statistically significant 
claim for either detection or the lack of the GZK feature. 

A comparison of the spectra obtained from AGASA and HiResI shows a sys- 
tematic shift of the two data sets, which may be interpreted as a systematic 
error in the relative energy determination of about 30%. If no correction for 
this systematic shift is carried out, the AGASA data are best fit by an injec- 
tion spectrum E' 1 with 7 = 2.6 for energies above 10 19 eV. The fit steepens 
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Fig. 8. Simulated spectra for the best fit injection spectrum with 7 = 2.6. Upper 
panel shows simulations for the AGASA event statistics after correcting for energy 
by an overall shift of -15%. The lower panel shows the fluctuations expected for the 
event statistics of HiResI after shifting HiResI energies by +15%. The shifted data 
for AGASA (grey circles) and HiResI (dark squares) are shown in both panels. 

to 7 = 2.8 when considering events down to 10 18 ' 5 eV. For HiResI, the best 
fit is between 7 = 2.7 and 7 = 2.8 for events with energy above 10 19 eV and 
7 = 2.6 for events above 10 18 ' 5 eV. With these best fits to the injection spec- 
trum the AGASA data depart from the prediction of a GZK feature by 2.6<r 
for 7 = 2.6. The HiRes data are fully compatible with the prediction of a GZK 
feature in the cosmic ray spectrum. The fit to the data with energy above 10 19 
eV is probably less susceptible to contamination by a possible galactic flux. In 
this case the AGASA departure from the expected GZK feature is, as stressed 
above, at the level of about 2.6a. Taking into account the uncertainties de- 
rived from the simulations, attributed to the discreteness of the photo-pion 
production and to cosmic variance, this discrepancy becomes even less signif- 
icant (~ 2.3a). It is clear that, if confirmed by future experiments with much 
larger statistics, the increase in flux relative to the GZK prediction hinted by 
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AG AS A would be of great interest. This may signal the presence of a new 
component at the highest energies that compensates for the expected sup- 
pression due to photo-pion production, or the effect of new physics in particle 
interactions (for instance the violation of Lorentz invariance or new neutrino 
interactions). 

Identifying the cause of the systematic energy and/or flux shift between the 
AG ASA and the HiRes spectra is crucial for understanding the nature of UHE- 
CRs. This discrepancy has stimulated a number of efforts to search for the 
source of these systematic errors including the construction of hybrid detec- 
tors, such as Auger, that utilize both ground arrays and fluorescence detectors. 
A possible overestimate of the AGASA energies by 15% and a corresponding 
underestimate of the HiRes energies by the same amount would in fact bring 
the two data sets in agreement in the region of energies below 10 20 eV. In 
this case both experiments are consistent with a GZK feature with large error 
bars. The AGASA excess is at the level of 1.7a (lAa if the observational un- 
certainties are also taken into account). Interestingly enough, the correction 
by 15% in the error determination implies that the best fit injection spectrum 
becomes basically the same for both experiments (E~ 2S ). 

With the low statistical significance of either the excess flux seen by AGASA or 
the discrepancies between AGASA and HiResI, it is inaccurate to claim either 
the detection of the GZK feature or the extension of the UHECR spectrum 
beyond Eqzk at this point in time. A new generation of experiments is needed 
to finally give a clear answer to this question. In Fig. 9 we report the simulated 
spectra that should be achieved in 3 years of operation of Auger (upper panel) 
and EUSO (lower panel). The error bars reflect the fluctuations expected in 
these high statistics experiments for the case of injection spectrum E~ 2£ . (Note 
that the energy threshold for detection by EUSO is not yet clear.) It is clear 
that the energy region where statistical fluctuations dominate the spectrum 
is moved to ~ 10 20 6 eV for Auger, allowing a clear identification of the GZK 
feature. The fluctuations dominated region stands beyond 10 21 eV for EUSO. 
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